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ABSTRACT 

Global, general relativistic magnetohydrodynamic (GRMHD) simulations of 
nonradiative, magnetized disks are widely used to model accreting black holes. 
We have performed a convergence study of GRMHD models computed with 
HARM3D. The models span a factor of 4 in linear resolution, from 96 x 96 x 64 to 
384 X 384 X 256. We consider three diagnostics of convergence: (1) dimensionless 
shell-averaged quantities such as plasma /3; (2) the azimuthal correlation length of 
fluid variables; and (3) synthetic spectra of the source including synchrotron emis- 
sion, absorption, and Compton scattering. Shell-averaged temperature is, except 
for the lowest resolution run, nearly independent of resolution; shell-averaged 
plasma /3 decreases steadily with resolution but shows signs of convergence. The 
azimuthal correlation lengths of density, internal energy, and temperature de- 
crease steadily with resolution but show signs of convergence. In contrast, the 
azimuthal correlation length of magnetic field decreases nearly linearly with grid 
size. We argue by analogy with local models, however, that convergence should be 
achieved with another factor of 2 in resolution. Synthetic spectra are, except for 
the lowest resolution run, nearly independent of resolution. The convergence be- 
havior is consistent with that of higher physical resolution local model ( "shearing 
box" ) calculations and with the recent nonrelativistic global convergence studies 
of Hawley et al. (2011). 
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Introduction 



The numerical study of black hole accretion flows has advanced significantly in the 
last decade. The advent of techniques for numerically solving the equations of general 
relativistic magnetohydrodynamics (GRMHD) has enabled self-coi isistent global modeling 



of accretion d riven by the magneto-rotational instability (MRI) ( iBalbus fc Hawleyl Il991 



Gammid 120041 ) o nto rotating black holes. Qualitative aspects of these simulations are code 
independent (e.g. iDe Villiers fc Hawleyl l2003l : iGammie et al.ll2003uAnninos et al.ll2005l ). but 
quantitative variations raise the question of numerical convergence. Recent work has shifted 
focus from dynamical properties of the accretion flow to simulated o bservations that can 



potentially constrain parameters for particular sour ces such as Sgr A* (IDolence et al.l 12009 



Moscibrodzka et al. 


2010 


(Shcherbakov et al. 


2010) 



Dexter et al.l l2009l |2010|) , including polarized radiative transfer 



to assess quantitative convergence of the underlying GRMHD simulations. 

Convergence studies of global accretion models are computationally expensive. An alter- 
native is to use a local model with shearing box boundary conditions to study the dynamics of 
MRI driven turbulence. These are simpler in the sense that there are fewer free parameters, 
and cheaper in that numerical resolution can be focused on a few correlation volumes ~ H^, 
where H is the disk scale height. The local model has for decades been a key theoretical tool 
for probing astrophysical disks (e.g. iGoldreich fc Lynden-Belll Il965l : iGoldreich fc Tremaine 



19781 : iNarayan et al.lll987l ) coupled to the s^ 



learing box boundary conditions has been widely 



used for the study of magnetized dis 


cs (e.g. 


Hawlev & Balbusll991. : 


L992; 1 


iawlev et al. 


1995, 


1996; 


Stone et al. 


1996; 


Sano et al. 


200^ 


1; 


Hirose et al. 20061: Fromane & Papaloizou 


20071; 


Fromane et al. 200?!: Guan et al. 20091: 1 


3avis et al.l 2010; Fromang 


2010 


Guan & Gammie 


2011; 


Simon et al. 


boil) 





Shearing box models have been integrated (1) with or without a mean magnetic field; 
(2) with or without stratification; (3) with or without explicit dissipation; (4) with and with- 
out explicit treatment of energy transport. There are now dozens of shearing box studies 
that treat aspects of this problem. The only models that c learly do not converge are un - 
stratified, zero-net field models without explicit dissipation (IFromang fc Papaloizoul 120071 ) . 
These models have a ma gnetic field correlation length th at decreases proportional to the grid 
scale ( Guan et al. 2009). But with explicit dissipation ( Lesur fc Longarettj 20071: Fromang 
2010h. a mean field Mawley et al. I I1995I : iGuan et al.ll2009l ) , or stratification (jDavis et al.lboiol: 



Simon et al.ll201ll ). the models do converge. One of the best resolved studies is iDavis et al. 
(120101 ). who convincingly demonstrate convergence of a stratified, isothermal, zero explicit 
dissipation model with the athena code at a physical resolution of up to 128 zones per 
scale height H. These stratified local models are physically closest to global simulations 
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(e.g. iHirose et al.ll2004l ). which are dominated by toroidal magnetic field. Local studies have 
shown, therefore, that with sufficient resolution numerical studies of MRI-driven turbulence 
can converge. 

Local models can focus on a few H^^ while global simulations must contain many H 
Are any of the dozen or so global disk models (e.g . 



Brandenburg 



1996 



Matsumoto et al. 



1996 



De Villiers fc Hawlev 2003: Gammie et al. 



2003; De Villiers et al. 2003:lMcKinnev fc Gammie 



Fragile et al.lboOTt 'Seckw ith et al.ll2008 



2004 



20081: 



Gammie et al.ll2004: 



Beckwith et al. 



Penna et al.ll2010l 



McKinnev 



2009 



2006 



Beckwith et al.ll2011 



Fragile et al]l2009l: IPragile fc Meieni2009i:lNoble et al. 



Flock et al.ll2011 



Hawley et al. 



2011 



ers) co n verged? And are syn t hetic observations based on gl obal models ( e .gJPexter fc Fragile 



( 2011 ): iHilburn et al.l fl2010f): iMoscibrodzka et al.l (120101 ): IPexter et al.l ( |2009[ ): iNoble et al. 



Shafee et al, 



2009 



2010 



and many oth- 



( 120071 ) : ISchnittman et al.l ( 120061 ): Dolence et al. 2011 in prep. ) sensitive to resolution? While 



some authors have inclu ded limited resolution studies (e.g. IShafee et al.ll2008l : iNoble et al. 



2010l : iPenna et al.l 120101 ). the answer is not yet clear. 



The first sys t emati c convergence test of a global black hole accretion simulation was done 
by iHawley et al.l (120111 hereafter HGK), using a zeus type code to simulate an HjR ^0.1 
disk in a pseudo- Newtonian potential. HKG simulate a 7r/2 wedge in azimuth, varying 
resolution around a fiducial 256 x 288 x 64 (r, z, cf) in cylindrical coordinate). After reviewing 
local model simulations and global nonrelativistic models HGK find that a minimum of 10 
cells per vertical characteristic MRI wavelength is required for convergence (HGK's Qz'-, e.g. 
Sano et al.l 120041 ). and 20 cells per azimuthal MRI wavelength (H GK's Qs)- They co nclude 
that most global simulations to date are far from resolved, except iNoble et al.l (120101 ) which 
used barely adequate poloidal resolution. 

In this paper we study the same convergence problem considered by HGK, but (1) in 
relativistic MHD and (2) using slightly different diagnostics. We ask what resolution is 
required for convergence (if convergence can be achieved), and how the global resolution 
requirements are related to local models. We are also particularly interested in whether 
resolution inffuences the spectra calculated from the models in the weakly radiative limit. 
This requires a fully relativistic simulation since in weakly radiative accretion ffows much of 
the emission arises from plasma near or even inside the innermost stable circular orbit (ISCO) 
of a spinning black hole. At these radii the relativistic models incorporate the dynamics of 
the plunging region and strong lensing effects on the radiation field. 

A third contrast with HGK is that we simulate a full 27r in azimuth rather than 7r/2. Our 
experience suggests that there is structure in the disk in all azimuthal fourier components, 
with the most power in the m = 1 component. Models with small azimuthal extent have 
reduced field strength and therefore require higher physical resolution by the HGK Q criteria. 
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We proceed as follows. §2 describes the code and initial and boundary conditions. §3 de- 
scribes convergence of radial profiles of non-dimensional variables. §4 describes convergence 
of azimuthal correlation lengths. §5 describes convergence of simulated spectra calculated 
with a Monte Carlo code. §6 gives a brief summary. 



Simulations 



Throughout the paper, we follow the standard notation of iMisner et al.l ( 119731 ) and set 
GM = c = 1. We consider a test fluid (no self-gravity) in the Kerr metric with dimensionless 
spin a* = 1 — ^ 0.94. The governing GRMHD equations express conservation of particle 
number 

(P«^);^ = , (1) 

and conservation of energy-momentum 



T^;, = 

together with the source-free Maxwell equations 

= 



(2) 



(3) 



where u^, p, T^'^, and *F'^'^ are the fluid's four velocity, rest mass density, GRMHD stress- 
energy tensor, and dual of the electromagnetic field tensor, respectively. The equation of 
state is 

P=il- 1)m (4) 

where 7 = 13/9, appropriate for a coUisionless plasma with relativistic electrons and non- 
relativistic protons. 



We evolve the GRMHD equations using the HARM3D code f lNoble et al.l l2009l . I2OO6 



Gammie et al. I I2003I ). HARM3D is a conservative high-resolution shock-capturing scheme demon- 
strated to have second order convergence in space and time for smooth flow s. The zone- 
centered magneti c field is updated with flux-interpolated constrained transport (iGammie et al. 



2OO3I : iTothl I2OOOI ) which preserves a particular numerical representation of V • B = 0. For 
this study, we use piecewise parabolic interpolation for both fluxes and EMFs. 



The numerical grid is uniform in modified Kerr-Schild coordinates Xi, 
dGammie et"aDl2003h : 



Xi = Inr 
9 = 11x2 + /isin(27rx2) 

X3 = 



X2, and X3 

(5) 
(6) 
(7) 
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where r, 6, and (f) are the Kerr-Schild radius, colatitude, and azimuth, respectively. We set 
h = 0.35 to concentrate the grid near the equatorial plane. The grid extends from below 
the horizon to r = 40, [O.OlTvr, O.QSSvr] in colatitude, and [0, 2it) in azimuth. HARM3D sets 
a "floor" for density and internal energy to avoid numerical problems that arise when those 
values are low: pmin = 10~V~^/^ and Umin = 10~^r~^/^. 



The initial condition is an equilibrium, prograde torus (IFishbone fc Moncrieilll976[ ) with 
inner edge at r = 6, pressure maximum at 12, and outer edge at 40. To make the torus 
unstable to MRI, it is seeded with weak poloidal magnetic field whose vector potential is 



A, 



C{p/P max 

0.2) i{A^>0 
if < 



1 



where C is a constant and Pmax is the maximum initial density. This gives dipole field line 
loops that run parallel to density contours. The field strength is normalized so that the ratio 
of the maximum gas pressure to maximum magnetic pressure (3 is 100. Small perturbations 
are introduced into the initial conditions to seed the MRI. The density and magnetic field 
lines are shown in Figured] for the initial conditions and for a later snapshot of the turbulent 
accretion flow. 

The models have outflow boundary conditions at the inner and outer radial (a;i) bound- 
aries and periodic boundary conditions in the azimuthal (xs) direction. The remaining (X2) 
boundaries are offset slightly from the pole, so the grid excludes a narrow cone around each 
pole. This avoids having the last polar zone control the timestep via the Courant condition 
because the polar zones become narrow in X3 (the computational expense is proportional 
to A^^ if poles are included!). While this treatment is essential for a convergence study, it 
is difficult to implement an appropriate boundary condition on the cone. We consider two 
different polar boundary conditions. 

The first, "hard" boundary is a solid reflective wall. We manually set the flux through 
the boundary to zero, and adjust the EMF in the flux-ct routine to make the cutout com- 
pletely opaque to the magnetic field, since the field vectors are modified in the routine after 
setting the boundary condition. This boundary condition produces an unphysical relativistic 
flow in the grids along the polar cone, so in addition we force the poloidal velocity in the 
zones along the boundary to be zero. 

The second, "soft" boundary also models a reflective wall. The variables in the ghost 
zones are all copied from the first physical zone. The X2 components of the velocity and 
magnetic field are inverted across the boundary (as usual for reflecting boundaries), but 
this only zeros fluxes on the boundary to within truncation error. This version of the polar 
boundary condition permits some leakage of magnetic flux through the polar boundaries. 
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but does not produce unphysical flows along the boundary. 

We ran a low resolution simulation with no polar cutout to evaluate both boundary 
conditions. The results suggest that the difference between the boundary conditions does 
affect the evolution of the high latitude "funnel" region. The soft boundary condition, in 
particular, causes a steady drop in the funnel region magnetic flux. On the other hand, all 
three cases (hard, soft, and no cutout) exhibit remarkably similar disk evolution. 

Table 1. List of Rtiks 
Resolution Duration ( ''""'jf'" ) Polar Boundar}' T}'po 



96 X 96 X 64 16,000 Soft 

128 X 128 X 96 12,000 Soft 

192 X 192 X 128 10,000 Soft 

384 X 384 X 256 6, 000 Soft 

96 X 96 X 64 16,000 Hard 

128 X 128 X 96 12,000 Hard 

192 X 192 X 128 10,000 Hard 



Our runs have numerical resolution {N^^, N^^, N^^) = (96, 96, 64), (144, 144, 96), (192, 
192, 128), and (384, 384, 256). The runs last until = 16, 000 for 96 x 96 x 64, 12, 000 for 
144 X 144 X 96, 10, 000 for 192 x 192 x 128, and 6, 000 for 384 x 384 x 256. Each resolution 
is run for both the soft and hard polar-boundary conditions except the highest resolution 
case which is run only for the soft-polar boundary due to numerical expense. A list of runs 
is shown in Table 1. The runs required ~ 10^(A''a;^/384)^(t//6, 000) cpu hours on TACC 
ranger. 

Each simulation's initial data contains noise inserted in each zone with a random number 
generator. This noise seeds the growth of instabilities in the torus. Each run will therefore 
differ in the details of the evolution, but over long enough periods one expects the differences 
to average away. Nevertheless, because our runs have finite duration, we expect some "cosmic 
variance," and this noise from run-to-run variations is present in every measurement we use 
to evaluate convergence. 

To evaluate run-to-run variation, we have repeated each of the Nxj_ = 96 and N^-^ = 
144 runs 3 times, and have used the variance of these runs to attach error bars to our 
measurements. We find that large run-to-run variations are caused by "events" that last 
a non-negligible fraction of the simulation time. For example, the lowest resolution runs 
sometimes gather a large mass of plasma near the ISCO, then dumps it suddenly into the 
black hole. We have also observed a bundle of magnetic field directed opposite to the field 
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in the funnel merge into the funnel, leading to a large fluctuation in the run with resolution 
144 X 144 X 96 and hard-polar-boundary. While the nature, frequency, and origin of these 
events is still unclear (we have only a handful of runs) it appears that run-to-run variation 
decreases at higher resolution. 



3. Radial profiles of non-dimensional variables 

We will compare poloidally, azimuthally, and time averaged radial profiles of the flow 
variables for the different resolution runs. We take a density-weighted average to focus on the 
accretion flow within ~ if of the equatorial plane. The explicit expression for the averaged 
radial profile F{xi) for a variable / is 

F{xi) = '\ (9) 
^2 — n 

where 

fit, xi) = V) ' (10) 

is the density weighted poloidally and azimuthally averaged radial profile of the variable / 
and g = g{x) is the determinant of the metric. For our case, ((2:2)1, (2:2)2) = (0.01, 0.99) and 

((X3)l,(x3)2) = (0,27r). 

We compare only non-dimensional variables since dimensional variables depend on the 
accretion rate, which decreases in time as the initial torus is accreted by the black hole. 
Our choice of the non-dimensional variables are scaled electron temperature 6e = kTe/rUe 
(= mppg/{2m,p) if Te = Tp), and /3 = Pg/pe = {T - l)u/{b^/2), where = h^'b^, 

b^ = ^(g>^^ + u^'u,)B'' (11) 

7 



= -n^*F'"' where = {-^-l/g'\ 0, 0, 0) , (12) 

7 is the Lorentz factor of the flow measured in the normal observer's frame, and F, k, rUp 
and rrie, and Tp and Te are the adiabatic index, Boltzmann constant, proton and electron 
mass, and proton and electron temperature, respectively. When calculating (3 we average Pg 
and pb separately using equation [10] and take the ratio of the averages. This prevents zones 
with near-zero magnetic energy from dominating the average. 

Figure [2] shows the radial profile of (3 and 9^ calculated using equation [9] for all the 
runs. All time averages run from t = 4000 to the end of the run; at t = 4000 the disk 
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at r < 10 is in a steady state for all runs except for the lowest resolution model, which 
shows a clear upward trend in /3 over the entire run. The lowest (96 x 96 x 64) and medium 
(144 X 144 X 96) resolution runs are averaged over 3 runs with different initial seeds to reduce 
run-to-run variation. The figure shows profiles for both the hard and the soft polar boundary 
conditions described in §2. 

Figure [3] shows (3 and 9e plotted against radial resolution A^^.^ for r = 2.04 (ISCO) 
and 8. The soft- and hard-polar boundary results are shown as solid black and red lines, 
respectively. Most quantities vary sharply from A^^.^ = 96 to 144 and then far less at higher 
resolution. For example, the soft polar boundary models have /3(ISC0) = (11.6, 7.3, 7.8, 6.6), 
and ^e(ISCO) = (31,47,48,57) at the four resolutions. 

Notice that at resolutions greater than 144 x 144 x 96 there are only small quantitative 
differences between the hard- and soft-polar boundary conditions, as seen in Figure [2] and [31 
We conclude that the effect of the polar boundary conditions on the main, equatorial flow 
is small for these dimensionless variables. 

What part of the variations at A^^^^ > 144 is real variation with resolution, and what 
part is run-to-run noise? The error bars in Figure [3] show standard deviation of the three 
runs performed for the lowest (96 x 96 x 64) and medium (144 x 144 x 96) data points with 
different initial seeds. Error bars are not available for the higher resolution data points due to 
computational expense. The size of the error bars is comparable to the differences between 
models run with different resolution. One might hope to gain additional information by 
measuring, e.g., (3 at several radii and averaging the trend with resolution, but, interestingly, 
the entire radial profile varies in a correlated way. Nevertheless Figures 2 and 3 show a clear 
trend of decreasing f3 and 6e with increasing resolution. It seems likely, therefore, that there 
is a genuine but weak trend with resolution. 



4. Correlation lengths 

We have looked at one-point statistics for non-dimensional variables. What about two- 
point statistics, which measure the spatial structure of the turbulence, and in particular the 
correlation length? The correlation length is a natural measure of the outer scale of the 
turbulence, and should be resolved and independent of resolution in a converged simulation. 

We consider only the azimuthal correlation length, as this is most straightforward to 
compute, and is most often under resolved in global simulations (HGK). The correlation 
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function at radius r on the equatorial plane is 

R{^) = r ^fiMSfi^o + <^)#o , (13) 

where 6f is deviation from average value of variable / at r. In practice, we average R in 
small area rArA6 across the equatorial plane, normalize, and average in time: 

R{r,(j),t)= R{r,e,(j),t)rdrde/{rArAe) (14) 

R(r,(f))=[ R{r,(j),t)/R{r,0,t)dt . (15) 
Jti 

Note that the correlation function for magnetic field is defined as 

i?(0) = / 56^(0o)%(0O + 0)^^00 , (16) 

where 6^ is defined in §2. Then 

A : i?(r,A) = i?(r,0)/e . (17) 

is the correlation length at radius r. 

Figure H] shows the azimuthal correlation length for density p, internal energy u, mag- 
netic field b, and 6e for all runs. Evidently the correlation lengths (angles) are nearly in- 
dependent of r, except close to the outer boundary where the models are not in a steady 
state. The correlation length varies between about 0.27r at the lowest resolution to O.lvr at 
the highest resolution for =dl variabte except i. Since H/r ~ 0.30 for all models over a wide 
range in radius (Figure E]), this corresponds (assuming flat space geometry) to 1 to 2 vertical 
scale heights. 

The non-dimensional resolution A/A0 ~ 12(A/(///r))(A^^^/384) where A0 = 2tx/Nx^, 
is marginal even for our highest resolution simulation. For 6, the correlation length of the 
highest resolution is smaller than that for any other variable. The magnetic field structure 
is underresolved. 

Figure [5] plots correlation length against resolution at the ISCO for the same variables 
as in Figure HI here red is the hard polar boundary and black is the soft polar boundary. 



^The scale height at each radius is defined as average of JJ^ {9 — 7r/2)^pd6'/ JJ^^ pd9 and /J^2 ° {9 — 
-K jlY pdB I J^i-f" pd9 where Oq is colatitude angle of the cutout = O.OITtt. 
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The dotted lines show how the correlation length would vary if it were fixed at 2, 5 and 10 
grid zones. 

For p, u, and 6^ (the nonmagnetic variables) the correlation length is ~ 5 grid zones 
for the two lowest resolution simulations. At higher resolution- N^-^ = 192 and 384- the 
correlation length increases to > 10 grid zones, and the slope of the change in correlation 
length with resolution decreases. This suggests that for the two highest resolution runs some 
structures in the turbulence are beginning to be resolved. 

For b, on the other hand, the correlation length decreases nearly proportional to the grid 
scale, with the correlation length fixed at around 5 grid zones per correlation length. There 
are small signs of an increase at the highest resolution, but in light of run-to-run variations 
the significance of this increase is marginal at best. The outer scale for the magnetic field is 
not resolved. 

For all variables the correlation lengths for hard and soft boundary polar conditions are 
consistent. Evidently the polar boundary does not infiuence the structure of turbulence in 
the equatorial disk. 

How do these correlation lengths correspond to those found in local model simulations? 
Guan et al.l (|2009[ ) found in their unstratified shearing box model that the three dimensional 
correlation function was a triaxial ellipsoid elongated in the azimuthal direction and tilted 
int o trailing ori e ntatio n. The relationship between our azimuthal correlation length A;, and 
the iGuan et al.l ( 120091 ) results is 

-1/2 



A 



cos^ ft 



tilt 



\2 

maj 



+ 



sin^ e 



tilt 



(18) 



where Quit ~ 15deg is the tilt angle of the correlation ellipse, and Xmaj, ^min are the major 
and min or axis of na a gnetic correlation lengths. For the best resolved net azimuthal field 
model in Guan et al. (120091 ) (y256b, which like our global models saturates at /9 ~ 20), this 
implies A ~ 0.17H ~ 0.05 rad, or O.OlGvrrad. Therefore, it is surprising that correlation 
length as large as ~ 0.3 rad ~ H are measured in our model for the nonmagnetic variables. 



Davis et al.l (120101 ) have computed correlation lengths in stratified, isothermal models 



with zero net flux. In a model run with athena at a resolution of 64 zones per scale height, 
the implied azimuthal correlation length (averaged over —H < z < H) fo r the magnetic 
field is slightl y larger than in the uns tratified models of I Guan et al.l (120091 ) . about 0.23H, 
or 0.027rrad. iGuan fc Gammid (l201ll ) have also run stratified, isothermal models at lower 
resolution with a zeus type code. They find an implied azimuthal midplane correlation 
length (similarly averaged) for the magnetic field that is even larger, about 0.9H, or 0.097r 
radians. Since correlation length decreases with increasing resolution it is possible that 
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Guan fc Gammid (120111 ) are not resolving the correlation lengt h, and that a t high er resolu- 
tion the correlation length would be closer to that measured by iDavis et al.l (|2010[ ). 



The correlation length of our highest resolution run spans 0.6(if/r) to OA{H/r) from 
ISCO to r ~ 10 where the corr esponding i s 7 an d 16, respectively. This is larger than the 
stratifi ed shearing box results of lDavis et al.l (120101) but smaller than t hat of I Guan fc Gammie 
( 1201 ll ). To resolve the correlation length found in lDayis et al.l (20101) we wou ld need another 
factor of 2 in linear resolution. Note that recently iBeckwith et al.l ( 1201 ll ) found in their 
global thin disk MHD simulation that azimuthal correlation length to be about 1.3{H/r) by 
averaging \z \ _ < H and 5 < r < 11. Th i s is la rger than our result but also falls between 



Davis et al.l (l2010f ) and lGuan fc Gammid ( 120111 ). 



5. Spectra 

An int eresting application of GRMHD models is to simulate observations of sources such 



as Sgr A* (iDolence et al. 



2009 



20091 . boioUPexter fc Fragilel 



Moscibrodzka et al. 2010; Hilburn et al. 2010; Dexter et al. 



201 ll ). Are the simulated spectra converged? 



The dynamical models underlying the spectral models are run with zero cooling, and 
the spectra are produced in a post-processing step. This is self-consistent as long as the 
fiows are advection dominated: the accretion timescale is much shorter than the cooling 
timescale. We calculate the er nergent radiation us ing grmonty, a general relativistic Monte 
Carlo radiative transfer code ( iDolence et al.ll2009l ). 



grmonty makes no symmetry assumptions and includes synchrotron emission, absorp- 
tion, and Compto n scattering. Using the rest-frame emissivity for a hot, thermal plasma 
( Leung et al.lboil" ) the code produces Monte Carlo samples of the emitted photons- "superphotons" 
that carry a "weight" representing the number of photons per superphoton. The superpho- 
tons follow geodesies, with weight varying continuously due to synchrotron absorption. They 
also Compton scatter and produce new, scattered superphotons with weight proportional to 
the scattering probability. We use a "fast light" approximation, where for each snapshot of 
simulation data a spectrum is created by treating the fiuid variables as if they were time- 
independent. This approximation is excellent for the time-averaged spectra we consider here. 
Superphotons that reach large radius are collected in poloidally and azimuthally distributed 
bins, and each bin p roduces a spectrum. A complete description of the code is given in 
Dolence et all toO^ . 



To compare runs we generate spectra for 200 — 1200 time slices (depending on the 
length of the run) and time-average them. The spectrum of each time slice is produced from 
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azimuthally averaged bins that extend from 0.127r < \9 — 7r/2| < O.lSvrrad with respect to 
the equatorial plane. 

We modify the simulation-provided data in one respect before calculating the spectrum. 
The quality of the non-magnetic fluid variable integration in the funnel region is poor due 
to truncation error. In particular the temperature can be high {9^ > 10^) and the particle 
density is determined entirely by a density floor in HARM3D. We therefore zero the emissivity 
if b'^/p > 1 to avoid contaminating the spectrum with possibly unphysical emission. 

It is necessary to fix a mass, length, and time unit to generate a radiative model. The 
combination GMbh sets a length and time scale but not a mass scale because the mass of 
the accretion flow is negligible in comparison to the black hole. We set Mbh = 4.5 x IO^Mq, 
comparable to the mass of SgrA*. The mass unit for the torus M. is still free; we set it 



so th at the 1.3 mm flux matches the observed flux from Sgr A* of ~ 3.4 Jy (iMarrone et al. 



20061). 



We want to model emission from a statistically stationary accretion flow. Because we 
start with a finite mass torus and it accretes over time, however, there is a steady decrease in 
density, field strength, accretion rate, etc., as the simulation progresses. We scale away this 
long term evolution using a smooth model, as follows. We set the mass unit Ai = Mos{t) 
where Mq is a constant and s{t) is a two-parameter scaling function. Then 



Punit = M/{ ) Uunit = PunitC B unit = W AlT punit , (19) 



C 



or expressing with s(t), 



Punit = Pos{t) Uunit = Uos{t) Bunit = Boy/s{t) (20) 

where they are the unit mass density, internal energy, and magnetic field strength, respec- 
tively, and po, Uo, and Bq are constants. Conversion from the simulation unit to the cgs unit 

IS, e.g. Pegs PsimPunit- 

The scaling function we employ has a form 

^ Ar5/3expf-^^ (21) 



s{t) "V t 

where A and are free parameters determined by a fit to the numerical evolution. The 
form comes from fitting 1-d relativistic viscous disk models (see Dolence et al. 2011 in prep, 
for more complete discussion). Notice that without this time-dependent scaling procedure, 
or with a different scaling procedure, the spectra would vary systematically over the course 
of the simulation. The spectra would also differ systematically with resolution because the 
plasma /3 varies with resolution. 
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We fit for A and tlie viscous timescale ti, from simulation data after a quasi-steady state 
has been readied, typically from t = 2000 onwards. A sample fit to M, for the 192 x 192 x 128 
run, is shown in Figure [71 The variance of the normalized accretion rate decreases with 
resolution, that is, at higher resolution the fluctuations are smaller and equation |2I] gives an 
increasingly good fit. The maximum of the normalized accretion rate is nearly independent of 
resolution, when models with different resolution are compared over the same time interval. 

Broadband, time-averaged synthetic spectra are shown in Figure [HI The mass unit of 
the torus is fixed by the condition that /i/(230 GHz) = 3.4 Jy for a Sgr A* model measured 
at the solar circle. The shape of the spectrum is broadly similar at all resolutions for both 
polar boundary conditions. 

Figure [HI shows flux density plotted against resolution in the infrared (3.8 fim) and X-ray 
(integrated from 2 keV to 8 keV) where most of the emission is from direct synchrotron and 
single Compton scatterings, respectively. Some of the variation is likely due to run-to-run 
variation, as indicated by the error bars on the N^-^ = 96 and N^^^ = 144 models. The flux 
varies with resolution by less than about 50% at infrared and 30% at X-ray for A^^,^ > 144. 
The spectra therefore appear remarkably consistent and independent of resolution, at least 
for the M and M appropriate to Sgr A*. 

In a sense this is not surprising, because (1) our normalization procedure removes much 
of the variation that might arise from the decrease of /3 with resolution, and (2) the tem- 
perature is very well converged. The combined effect of the fixed flux normalization and 
the variation with resolution is to strengthen the magnetic field slightly and move the syn- 
chrotron peak slightly further into the infrared. This is echoed in the first Compton bump 
in the X-ray, which is forced to slightly higher energy by the increase in infrared input pho- 
tons. While we have demonstrated this for only a single set of the model parameters (M, 
/i,(230 GHz)), exploration of slightly different models with similarly consistent results shows 
that this is not a unique case. 

6. Summary 

We have investigated convergence of global GRMHD simulations of hot accretion flows 
onto a black hole and the emergent spectrum. We have run GRMHD simulations for four 
different resolutions, 96 x 96 x 64, 144 x 144 x 96, 192 x 192 x 128, 384 x 384 x 256 in 
spherical-polar coordinates. We have probed convergence using three diagnostics: time- 
averaged radial profiles of nondimensional quantities (plasma /3 and electron temperature 
6'e); azimuthal correlation lengths for several variables including the magnetic field; and 
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artificial spectra generated with a Monte Carlo code. 

For most of our diagnostics there are substantial differences between the lowest (96 x 
96 X 64) and next lowest (144 x 144 x 96) resolution, and relative minor changes at higher 
resolution. Run-to-run variations in the lower resolution models tend to be larger than the 
differences between the higher resolution (192 x 192 x 128 and 384 x 384 x 256) models. 

We find that the magnetic correlation length is not converged. It decreases nearly lin- 
early with resolution, with the number of grid cells per magnetic correlation length fixed 
at ~ 5, although we do see a slight increase as resolution increases. Comparison with local 
model/shearing box simulations suggests that the turbulence does not change quahtatively 
at higher resolution. Such comparisons also suggest that another factor of ~ 2 in linear reso- 
lution (costing about 1.6 x 10^ cpu-hours) would resolve the azimuthal magnetic correlation 
length. None of the existing simulations (local or global) resolve scales more than a factor 
of ~ 4 smaller than the correlation length (particularly the minor axis correlation length, 
which is oriented nearly along the radial unit vector and which we have not investigated 
here). If we identify the correlation length with the outer scale of MRI driven turbulence, 
as seems reasonable, then none of these models have a resolved inertial range. 

On the other hand, time-averaged synthetic spectra based on the GRMHD models, 
with parameters fixed to match Sgr A*, are remarkably reproducible from resolution to 
resolution. This suggests that simulated observations from existing simulations have some 
predictive power. We think it likely that the leading source of error in the high resolution 
radiative models is now related to the underlying physical model (particularly the fluid model 
treatment of the plasma, and the absence of conduction) rather than the finite resolution of 
the models. 

A similar convergence study has been conducted by HGK for nonrelativistic global 
models. It is worth asking whether our models are converged according to the dimensionless 
resolution Q, the ratio of most unstable MRI wavelength!^ to the grid cell size in the azimuthal 
and vertical direction. In the azimuthal direction. 



{Qy or Qif, in HGK's notation), where Cg ~ HQ is the sound speed. This gives ^ 22 and 




(22) 



(23) 



^Although Q is well defined, the background state is turbulent and there are no well defined linear MRI 
modes. 
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> 10 for Nx-^ = 384 and 192, respectively, for all radii less than 10. In the vertical direction. 



{Qz in HGK's notation) where A6 is the zone size in Kerr-Schild coordinates at the midplane. 
Since I^s/Esl is usually ~ 3 - 10, this gives Qa > 70 and > 30 for A^^^, = 384 and 192, 
respectively, for all r < 10. The required Q values to resolve the characteristic wavelength 
ci^re Q2 ^ 20 — 60 and Qs ^ 6. Hence, MRI in the toroidal direction is resolved but not in 
the poloidal direction in these runs according to HGK's Q criterion. 

To summarize our findings in the form of guidance for future simulators: (1) the reso- 
lution 96 X 96 X 64 is too low. The convergence measurements differ by factors of several 
from the highest resolution runs, and the magnetic field weakens steadily in a relative sense 
(/3 increases) over the course of the run; (2) the resolution 144 x 144 x 96 shows early signs 
of convergence except for the correlation length of the magnetic field; (3) the resolution 
192 X 192 X 128 and 384 x 384 x 256 differ relatively little from each other and show signs of 
convergence in the azimuthal correlation lengths, the temperature, and spectra, but not in 
the correlation length of magnetic field; (4) the observed trends with increasing resolution 
(to the extent that they are significant at the highest resolution) are that (3 decreases, 9e in- 
creases, correlation lengths decreases, and IR and X-ray fluxes increase relative to millimeter 
fluxes, which we use to normalize the spectrum. 
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Fig. 1. — Poloidal slices of the initial and turbulent state of the global simulation. The 
pseudo-color is showing scaled logarithmic density and black lines are the initial magnetic 
field lines. 
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Fig. 2. — Radial profile of plasma f3 (upper row) and electron temperature 6e (lower row) for 
each resolution. The columns are for the soft-polar-boundary (left) and hard-polar-boundary 
(right). 




Fig. 3. — Plasma /3 (left) and electron temperature 6^ (right) plotted as a function of res- 
olution at the ISCO (r = 2.04) and r = 8. The black lines are for the soft-polar-boundary 
and the red lines are for the hard-polar-boundary. 
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Fig. 4. — Azimuthal correlation length as a function of radius for each resolution. Prom 

the top panel, density (p), internal energy (m), magnetic field (&), and electron temperature 
(6'e). The left column is for the soft-polar-boundary and right column is for the hard-polar- 
boundary. 
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Fig. 5. — At ISCO, azimuthal correlation length of density (A^), internal energy {\u), mag- 
netic field (Afo), and electron temperature {9e) are plotted as a function of resolution. The 
black lines are for the soft-polar-boundary and the red lines are for the hard-polar-boundary. 
Black dotted lines show a correlation length of 2, 5, and 10 grid cells, to which correlation 
length size of 2, 5, and 10 grids correspond at each resolution in azimuthal direction. 
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Fig. 6. — Radial profile of the scale height H/r for the runs with the soft polar boundary. 
The runs with the hard polar boundary have similar profiles. 
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Fig. 7. — Time evolution of accretion rate for the run 192 x 192 x 128 with hard polar 
boundary. Dotted line is the actual accretion rate and the solid line is a fit of the form 
shown in equation fl2Tl) . 
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Fig. 8. — Spectra for each resolution. Flux is fixed to 3.4 Jy at 1.3 mm shown by the vertical 
solid line. The left plot is for the soft-polar-boundary and right plot is for the hard-polar- 
boundary. 




Fig. 9. — Infrared flux density (3.8 fim, left) and X-ray flux (integrated from 2keV to 8keV, 
right) as a function of resolution. The black lines are for the soft-polar-boundary and the 
red lines are for the hard-polar-boundary. 



